Three-dimensional simulations of distorted Black Holes. I. Comparison with 

axisymmetric results. 
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We consider the numerical evolution of black hole initial data sets, consisting of single black holes 
distorted by strong gravitational waves, with a full 3D, nonlinear evolution code. These data sets 
mimic the late stages of coalescing black holes. We compare various aspects of the evolution of ax- 
isymmetric initial data sets, obtained with this 3D code, to results obtained from a well established 
axisymmetric code. In both codes we examine and compare the behavior of metric functions, appar- 
ent horizon properties, and waveforms, and show that these dynamic black holes can be accurately 
evolved in 3D. In particular we show that with present computational resources and techniques, 
the process of excitation and ringdown of the black hole can be evolved, and one can now extract 
accurately the gravitational waves emitted from the 3D Cartesian metric functions, even when they 
carry away only a small fraction (<< 1%) of the rest mass energy of the system. Waveforms for 
both the I = 2 and the much more difficult I = 4 and I — 6 modes are computed and compared with 
axisymmetric calculations. In addition to exploring the physics of distorted black hole data sets, 
and showing the extent to which the waves can be accurately extracted, these results also provide 
important testbeds for all fully nonlinear numerical codes designed to evolve black hole spacetimes in 
3D, whether they use singularity avoiding slicings, apparent horizon boundary conditions, or other 
evolution methods. 
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I. INTRODUCTION 

In a recent paper we showed that full 3D, Carte- 
sian coordinate simulations of (axisymmetric) distorted, 
dynamic black holes can be performed, and through a 
comparison with results obtained with an axisymmetric 
code, that the waveforms can be accurately obtained. 
Although with present resolutions the length of the sim- 
ulations is limited to about t w AQMadm (due to the 
use of singularity avoiding time slicings, as detailed be- 
low), these results represent a step forward towards ob- 
taining useful waveforms from 3D numerical calculations 
that will be needed for detecting and analyzing signals 
buried in the data collected by gravitational wave obser- 
vatories. The need for an increasingly complex sequence 
of such simulations leading to realistic coalescence calcu- 
lations is great, as gravitational wave observatories may 
detect such waves in about five years As black hole 
collisions are presently considered a most likely source of 
signals to be detected initially by these observatories, it 
is crucial to have a detailed theoretical understanding of 
the coalescence process that can only be achieved through 
numerical simulation In particular, it is most im- 

portant to be able to simulate accurately the excitation of 
the coalescing black holes, to follow the waves generated 
in the process, and to extract gravitational waveforms 



expected to be seen by detectors. 

These are very difficult calculations, as one must simul- 
taneously deal with singularities inside the black holes, 
follow the highly nonlinear regime in the coalescence pro- 
cess taking place near the horizons, and also calculate 
the linear regime in the radiation zone where the waves 
represent a very small perturbation on the background 
spacetime metric. The energy carried by these waves is 
typically found to be on the order of 10~ 6 — 1Q~ 2 Madm- 
At such low amplitude, both the generation and propa- 
gation of these signals are susceptible to small numerical 
errors inherent in numerical simulations. Furthermore, 
one must numerically extract these waves from the met- 
ric functions actually being evolved, which in 3D are 
usually the Cartesian metric functions (e.g., g xx , etc.). 
As these metric functions do not correspond directly to 
physical degrees of freedom, but rather to the gauge in 
which they are evolved, the determination of waves is 
provided by a complicated numerical procedure to iso- 
late gauge- invariant waveforms (described below), which 
can introduce further numerical errors. 

In axisymmetry, in the case of stellar collapse |B| or 
rotating collisionless matter || the evolution has been 
followed through black hole formation and emitted £ = 2 
waveforms have been studied. In this paper we focus on 
vacuum black holes, which exist in the initial data. Such 
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data sets are likely to be used in the fully general study 
of 3D black hole coalescence. Evolutions of vacuum black 
hole data sets have been carried out successfully for dis- 
torted black holes with rotation |?],|| and without ||,|l0| , 
and for equal mass colliding black holes , but with 

difficulty. These 2D evolutions can be carried out to 
roughly t — IOOMadm, where Madm is the ADM mass 
of the spacetime, although beyond this point large gradi- 
ents related to singularity avoiding slicings usually cause 
the codes to become very inaccurate and crash. In full 
3D, the simulations are more difficult as the numerical 
problems are exacerbated jl3| . We will show examples of 
this difficulty in sections below. 

These problems should be reduced through the use of 
apparent horizon boundary conditions (AHBC), that al- 
low one to avoid the singularities not by pathological slic- 
ing conditions, but by excising the region inside the hori- 
zon. To date, impressive successes have been achieved 
in ID Jl4-jr7j], and also with spherical black holes in 
3D [ p^|Jl^Jl9j . But these techniques have not yet been 
applied to dynamic, nonspherical black hole systems. A 
completely different approach, using characteristic evo- 
lution, has recently been developed very successfully for 
3D single black hole evolution, including distorted and 
rotating black holes pc[ , allowing extremely long term 
3D evolutions for the first time. But for truly dynamic 
black holes evolved in 3D using any techniques (Cauchy 
evolution, Characteristic evolution, or AHBC), careful 
extraction of waveforms has been only recently been in- 
vestigated in detail in Ref. Those results showed for 
the first time that the black hole ringdown dynamics, 
and the accompanying emission of low amplitude waves 
could be accurately computed in a full 3D simulation. 
Such results should be a useful testbed for techniques 
such as characteristic evolution and AHBC that enable 
long term black hole evolution. 

In this paper we follow up on the results of Ref. jjj] 
with a much more extensive and detailed examination 
of the evolution of axisymmetric black holes with a fully 
3D code. These black hole initial data sets correspond to 
Schwarzschild initial data with a nonlinear gravitational 
wave superimposed, creating a distorted Schwarzschild 
system that will oscillate and settle down to a spheri- 
cal hole with waves traveling out far from the hole. We 
show that with current techniques and computational 
resources available to 3D numerical relativity, distorted 
black holes can be evolved through the initial relaxation 
and part of the final ringdown period. We compare the 
evolution of metric functions and horizons obtained with 
an axisymmetric code written in spherical-polar coordi- 
nates using maximal slicing and a particular shift condi- 
tion, to full 3D evolutions with a very general 3D code 
written in Cartesian coordinates, capable of using differ- 
ent slicing and shift conditions. We show that the grav- 
itational waveforms can be followed and accurately ex- 
tracted from the numerical evolutions, even though they 



represent a small perturbation on the background space- 
time which is also being evolved. We also study the ex- 
traction of fully nonaxisymmetric modes with the present 
code. In principle, of course, with axisymmetric initial 
data one should not excite nonaxisymmetric modes, but 
the extent to which they may be excited due to numeri- 
cal effects in present codes needs to be investigated. We 
show below that numerical effects may excite such modes 
to some degree. This will provide important information 
when evolving true 3D, nonaxisymmetric black hole data 
sets, so that true signals can be clearly distinguished from 
numerical artifacts. 

The extension of this work to full 3D black hole ini- 
tial data, for which no testbeds exist at present, is in 
progress and will be published in the next paper in this 
series. In particular, nonaxisymmetric modes, such as 
the I — 2, m — 2 mode expected to be important in the 
ringing radiation for rotating black holes at late times || , 
and therefore an important signal for gravitational wave 
observations, can now be studied |^l|-^3||. 

The plan of the present paper is as follows. In sec- 
tion U we briefly review the initial data construction. 
In section [II, we describe both the 2D and the 3D nu- 



merical codes used in the present simulations, and we 
study the numerical evolution of several representative 
distorted black hole data sets, discussing the behavior of 
metric functions, horizons, and waveforms extracted. In 
section [V we summarize and describe future work in this 



scries. 



II. INITIAL DATA 

In this section we describe the axisymmetric initial 
data describing the distorted black holes to be evolved. 



A. Mathematical construction and numerical 
implementation 

There is by now a large body of literature on black hole 
initial data for numerical relativity. For multiple black 
holes, the early axisymmetric data sets of Misner |Q] and 
Brill and Lindquist pq] were generalized to include spin 
and momentum by York and coworker s |26[|, culminating 
in the definitive work in 3D of Cook |27H29|1. An alter- 
nate and simpler technique for computing multiple black 
hole initial data was recently developed by Brandt and 
Briigmann |p0| . 

In this section we briefly describe the mathematical 
construction of initial data evolved in this paper, which 
constitute yet another family of black hole data sets for 
numerical relativity. These data sets are single, time 
symmetric, distorted black holes of the type studied in 
Refs. [ pl|j3^ ] as axisymmetric initial data, and their evo- 
lutions have been studied extensively in Ref. 0. These 
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data sets have also been extended to include rotation 
|0,||3^], but we will not consider rotating black holes 
here. We will make extensive use of these data sets in 
our studies presented in this paper. They consist of a 
single black hole that has been distorted by the presence 
of an adjustable torus of nonlinear gravitational waves, 
of a form originally considered by Brill p4| . which sur- 
round it. This is discussed in much more detail in Ref. 
pOf , which focuses only on the initial data, where we 
study a broad family of full 3D distorted black holes, 
with and without rotation. We restrict ourselves here 
to nonrotating and axisymmetric cases, as studied previ- 
ously in axisymmetry by [|TJ. The amplitude and shape 
of the torus of waves can be specified by hand, as de- 
scribed below, and can create very highly distorted black 
holes. For our purposes, we consider them as convenient 
initial data that create a distorted black hole that mim- 
ics the merger, just after coalescence, of two black holes 
colliding in axisymmetry |fl2|| . 

Following fjlfl , we write the 3-metric in the form orig- 
inally used by Brill |3: 



df 2 



^ (e 2 i (drj 2 + d9 2 ) + : 



(1) 



where rj is a radial coordinate related to the standard 
Schwarzschild isotropic radius, which we refer to here as 
r, by 



M 

r = — e '. 
2 



(2) 



Then the rj coordinate is related to the Cartesian coor- 
dinates by 



\J x 1 + y 2 + z 2 = e v . 



(3) 



(Here, and in what follows, we set the scale parameter 
M in |3l[] to be 2.) In spherical symmetry, if the "Brill 
wave" function q vanishes, the Schwarzschild initial data 
set results, as shown below, so for small values of the 
Brill wave one may regard this system as a perturbed 
black hole. For large waves, the black holes can become 
extremely distorted. 

Given a choice for the "Brill wave" function q, the 
Hamiltonian constraint leads to an elliptic equation for 
the conformal factor -0. (As these data sets are time 
symmetric, the momentum constraints are trivially sat- 
isfied. In Ref. (3^] we consider the non-time symmetric 
case.) The function q represents the gravitational wave 
surrounding the black hole, and is chosen to be 



(4) 



Thus, an initial data set is characterized by the parame- 
ters (a, b, w, n, c), where, roughly speaking, a is the am- 
plitude of the Brill wave, b is its radial location, w its 



width, and c and the even integer n control its angu- 
lar structure. Note that we have generalized the original 
axisymmetric construction to full 3D by the addition of 
the parameter c, but in this paper we restrict ourselves 
to c = for comparison with axisymmetric results. A 
study of full 3D initial data and their evolutions will be 
published elsewhere pjlpKj ]. 

If the amplitude a vanishes, the undistorted 
Schwarzschild solution results, leading to 



tp = 2 cosh ( I 



(5) 



We note that just as the Schwarzschild geometry has an 
isometry that leaves the metric unchanged under the op- 
eration rj — > —rj, our data sets also have this property, 
even in the presence of the Brill wave. The radial point 
rj = is called the throat, as it is the surface that con- 
nects the identical geometries of two universes connected 
by an Einstein-Rosen bridge, or wormhole. As discussed 
in ||l3| , [35|| , this isometry condition can also be applied 
during the evolution and in Cartesian coordinates as well. 
We will make use of this condition in the work described 
in this paper. 

As discussed in Ref. (30), we compute the initial data 
sets on a 3D grid in spherical coordinates, because the in- 
ner and outer boundaries are easier to handle. The outer 
boundary is a Robin condition, and the inner boundary 
is provided by the isometry condition described above. 
Once we have the conformal factor, we can specify the 
3-metric, the form of which is given in Eq. (g). The ex- 
trinsic curvature vanishes, as these are time-symmetric 
initial data sets. Because of the higher resolution of the 
spherical grid in the inner regions, we found that com- 
puting derivatives of the conformal factor on this grid 
and storing the interpolated values on the Cartesian grid 
leads to more accurate evolution. Thus, we compute the 
first and second derivatives of the conformal factor tp 
with respect to the spherical coordinates (77, 9, <p) on the 
spherical grid. 

For convenience, these manipulations are carried out 
on the spherical grid, but we perform evolutions in Carte- 
sian coordinates. We have found it adequate to inter- 
polate these quantities from the very fine 2D spherical 
grid onto the Cartesian grid using linear interpolation. 
With this information, we can perform simple coordinate 
transformations to convert these quantities to Cartesian 
coordinates obtaining: 
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using the coordinate transformation implicit in Eq. (^). 
Finally, ip and its Cartesian derivatives, computed via 
the coordinate transformation, are stored for use in com- 
puting "conformal" numerical derivatives, as described 
in Ref. Q. 
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Problem 


a 


a 


Madm/M 


Mah/M 


MRL 


1 


0.5 


2 


0.92 


0.86 


7.0 x W-'- 1 


2 


0.5 


4 


0.97 


0.92 


6.0 x 10~ 2 


3 


0.1 


4 


0.97 


0.96 


4.8 x 10~ 3 



0.50 



TABLE I. We give the parameters for the three initial data 
sets from whose evolution we will present results in this pa- 
per. The complete set of black hole initial data parameters 
are (a,b,w,n, c). The amplitude a and angular index n are 
given, while for all cases the Brill wave location and width pa- 
rameters are b = 0, w = 1, and c = indicating axisymmetric 
initial data. The quantity MRL is the maximum radiation 
loss possible during evolution for each data set, as discussed 
in the text. M is the scale parameter described in the text, 
taken to be 2 in this work. 



B. Characteristics of sample initial data sets 

In this section, we single out three initial data sets 
whose evolution and waveform emission are studied in 
detail in the following sections, and we present some of 
their physical characteristics before going on to evolu- 
tions. 

In Table |, for each initial data set we give the ADM 
mass Madm, the apparent horizon mass Mah, and the 
maximum radiation loss MRL, defined as in Ref. Bill as 



MRL 



Madm - M AH (t = 0) 
Madm 



(7) 



The MRL is an upper limit on the amount of energy 
the spacetime can emit as computed by considering the 
second law of black hole dynamics. The three cases stud- 
ied range from slightly distorted, with an MRL of about 
5 x 10~ 3 Madm, to highly distorted, with an MRL of 
7 x 10 -2 Madm ■ Problems 1 and 2 have the same distor- 
tion amplitude parameter, but different angular index n, 
leading to quite different ratios of £ = 4 to I = 2 content 
in the initial data. 

To give a visual picture of how distorted the black holes 
in these initial data sets are, we show the embedding of 
their apparent horizons. In Figs. [l], pi and [j] we show 
the apparent horizon embeddings for Problems 1,2, and 
3, respectively. As shown previously in Ref. [MJ, black 
holes of this type with positive distortion amplitude a are 
have prolate horizons. Problem 3 is slightly perturbed, 
and nearly spherical, while Problems 1 and 2 are notice- 
ably prolate. As in previous papers, we can character- 
ize the shape of axisymmetric black holes with a single 
parameter: the polar to equatorial circumference ratio 
C r . The horizons for Problems 1 and 2 are quite prolate, 
with C r = 2.7 and 2.4, respectively. Problem 2, with its 
angular index n = 4, actually has a "waist" , with nega- 
tive Gaussian curvature in the region near the equator. 
Axisymmetric (2D) evolutions of similar data sets were 
studied extensively in Refs. |^,|l^,^l|, and dynamics of 
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FIG. 1. We show the embedding of the apparent horizon 
for Problem 1, (a,b,w,n,c) = (0.5,0,1,2,0). This black hole 
horizon is quite prolate, with a polar to equatorial circumfer- 
ence ratio C r = 2.7 



the apparent and event horizons horizons of such data 
sets were detailed in 36-39 



These data sets will all go through oscillations, even- 
tually settling down to become spherical Schwarzschild 
black holes as energy is radiated away. Note that Prob- 
lems 1-3 form a sequence from most to least distorted, 
based on their horizon geometries, and also from largest 
to least possible radiation output, measured by their 
MRL. While Problem 3 is slightly distorted and a good 
candidate for study via perturbation theory, Problems 1 
and 2 are fairly strongly distorted and are not well de- 
scribed by a first order perturbative treatment [ p2"|]23]| . In 
the next section will study the full 3D evolution of these 
and other data sets in Cartesian coordinates. 



III. EVOLUTION OF AXISYMMETRIC INITIAL 
DATA SETS 

We now turn to the 3D evolution of the series of ax- 
isymmetric initial data sets described above. Compar- 
isons of these results with those from well established 
2D codes will help us confirm various aspects of the 3D 
results. 



A. Numerical codes 

1. 2D code 

The code used to produce the 2D results in this paper 
is based on the one described in p5[ . This code uses a 
logarithmic radial coordinate r\ related to the standard 
Schwarzschild isotropic radius r by r\ — In (2r/M), where 
M is a scale parameter. As all data sets evolved with 



4 



0.50 




0.20 : 



0.10: t 

0.00 F : 

0.00 0.10 0.20 0.30 0.40 0.50 

p/Va 

r/ AH 

FIG. 2. We show the embedding of the apparent horizon for 
Problem 2, (a, b, w, n, c) = (0.5, 0, 1, 4, 0). This horizon has a 
waist at the equator, with a polar to equatorial circumference 
ratio C r = 2.4 
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FIG. 3. We show the embedding of the apparent horizon 
for Problem 3, (a,b,w,n,c) = (0.1,0,1,4,0). This sightly 
prolate horizon is the least distorted of the three Problems 
studied. 



this code have equatorial-plane symmetry, the polar an- 
gle 9 was restricted to lie between and ir/2 to decrease 
computational time. All 2D results presented were ob- 
tained using 200 radial zones and a maximum radius of 
Vmax = 6, resulting in a radial grid spacing of Aij = 0.03. 
The number of angular zones, chosen to make AO « A77 
was approximately 50. 

The use of the logarithmic radial coordinate has the 
advantage of providing fine resolution near the throat of 
the black hole and also near the peak that develops in 
the radial metric function, while also allowing the outer 
boundary to be placed far from the hole. Disadvantages 
of this coordinate are that (i) The throat region remains 
extremely well resolved, even after the horizon has moved 
significantly away from the hole. Therefore much com- 
putational effort is wasted well inside the horizon where 
the lapse is typically near zero and the region is causally 
disconnected from the outside, (ii) The grid becomes 
very coarse outside the horizon in the radiation zone, 
because equal spacing in the rj coordinate leads to larger 
and larger spacing in the more physical r coordinate. Un- 
der these conditions, waves may be reflected back toward 
the black hole as they are scattered off of the coarse grid 
at larger radii, as discussed in ||J3^]. In fact, the effective 
Ar at the location of our wave extraction in the 2D code 
is approximately 0.45Af, which is typically much larger 
than the resolution of the 3D grid. 

Except for the tests performed with geodesic slicing, 
the slicing used in the 2D code was maximal slicing. In 
order to stabilize the code along the axis, the shift vector 
is chosen to maintain the metric function 7^0 = 0, thus 
keeping constant-coordinate lines perpendicular (7,^ and 
700 also vanish). These details are extensively discussed 
in Refs. 

2. 3D code 

a. The code. We developed a 3D code to study 
black holes and gravitational waves in Cartesian coor- 
dinates. This code (known as the "G" code) was applied 
to Schwarzschild black holes, where we showed that using 
singularity avoiding time slicings, a spherical black hole 
could be evolved accurately to t = 30 — 50M, depending 
on the resolution, location of the outer boundary, and 
the slicing conditions Beyond that time, the code 

generally crashes due to the unbounded growth of met- 
ric functions generated by singularity avoiding slicings. 
However, the focus of [jl3| was on spherical black holes, 
so no studies were made of black hole oscillations and 
the waves that would be generated in the process. It was 
shown that with spherical initial data some nonsphcri- 
cal behavior could be introduced by the Cartesian mesh 
and boundary conditions, the numerics of which could 
in principle generate spurious gravitational waves. This 
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3D code was then applied to the collision of two axisym- 
metric black holes (Misner data) ^(J, where we showed 
by comparison to 2D results that one could accurately 
track the merging of the horizons, and that the radiation 
emitted was qualitatively the same, but at that time the 
waveforms were not studied extensively. Building on the 
work presented in this paper and in Refs. p],pl[, a more 
detailed study of the Misner data in 3D, including the 
waveforms, is in preparation for publication elsewhere. 

The same code was simultaneously applied to the prob- 
lem of pure gravitational waves where many sys- 
tems were studied, from pure linear quadrupole waves 
to nonlinear waves, and their propagation on a Carte- 
sian mesh was studied. In that study it was shown that 
waves can be accurately evolved, although certain prob- 
lems with gauge modes in the "near linear" regime that 
can confuse the results were identified, along with strate- 
gies to deal with them. Finally, the code has also been 
used to study coordinate conditions in 3D numerical rel- 
ativity [||. 

Based on this "G" code, we have developed a new ver- 
sion using the same numerical techniques, having similar 
convergence properties, etc., but rewritten to take advan- 
tage of newer parallel computer architectures, notably 
the SGI/Cray Origin 2000. We use the shared mem- 
ory paradigm for parallelism in this code, avoiding the 
need to use a message passing library. With this model, 
we find both efficient single processor use and reason- 
ably good scaling beyond 64 processors on the Origin 
2000. The numerical techniques are as described in detail 
in We evolve the 3-metric and extrinsic curvature 

components with an explicit, staggered leapfrog scheme. 
For the algebraic slicing used to produce the results in 
this paper, described in the next section, the lapse func- 
tion is evolved in this scheme as well. 

b. Gauge conditions The 3D code has been written 
in a general way to accept general shift and slicing condi- 
tions. Although various shift conditions have been tested 
(e.g. minimal distortion shift as reported in Ref. |f43|j or 
an apparent horizon shift condition reported in Ref. [jl 3[j ) , 
in this work we restrict ourselves to the zero shift case. 
The application of different shift conditions to black hole 
spacetimes in 3D will be presented elsewhere. 

We have experimented with a number of different slic- 
ing conditions. In Refs. [l3pll| we discussed the use of 
geodesic slicing, maximal slicing and a host of algebraic 
slicing conditions for single and two black hole space- 
times. Geodesic slicing is useful only for performing cer- 
tain tests of a code on a black hole spacetime. In brief, 
the slicing condition that worked the best for the analytic 
Schwarzschild and Misner two black hole spacetimes was 
the evolution of an initially maximal lapse, which van- 
ishes on the throat of the black hole, with the "1+log" 
algebraic slicing. With this slicing condition we were able 
to produce accurate and stable evolutions for a period of 
time long enough to study the ringdown of excited black 



holes, and it is much more computationally efficient than 
maximal slicing. 

However, for the numerically generated distorted black 
hole initial data sets there is no appropriate analytic an- 
tisymmetric (i.e. vanishing on the throat) initial lapse to 
use, so we compute an antisymmetric maximal lapse nu- 
merically. Because it is inconvenient to compute and an- 
tisymmetric lapse across the black hole throat (which is a 
coordinate sphere) in the Cartesian coordinates with our 
elliptic solver, we compute the antisymmetric maximal 
lapse on the spherical grid and interpolate that onto the 
Cartesian grid. We have shown that this procedure when 
applied to the Schwarzschild case produces the same re- 
sults as when the analytic Schwarzschild lapse is used, 
although it is necessary to use at least quadratic inter- 
polation. This initial lapse is then evolved with the alge- 
braic "1+log" lapse as described in Ref. jl3|]. This is the 
slicing condition used in all runs in this paper. 

c. Computational domain and boundaries For the 
outer boundary, we presently hold the metric functions 
fixed in time. If the boundary is placed too close, this can 
have a serious effect on the behavior of waveforms and 
metric functions, and hence the largest possible compu- 
tational domain is desirable. For evolutions of 30 — 50M, 
typically we need to place the outer boundary at least as 
far as 20 — 25M. The inner boundary (on the throat) is 
provided by an isometry condition, as described above. 

Although the 3D evolution code is written without 
making use of any symmetry assumptions, the initial 
data we evolve in this paper have both equatorial plane 
symmetry and axisymmetry. Hence we save on the mem- 
ory and computation required by evolving only one oc- 
tant of the system. As shown in |L3j, this has no ef- 
fect on the simulations except to reduce the computa- 
tional requirements by a factor of eight. Even with such 
computational savings, these can be extravagant calcula- 
tions. Smaller simulations presented in this paper have 
resolutions of typically 100 3 , and can easily be computed 
in a few hours on a 32 processor Origin 2000 at AEI. 
The highest resolution results presented in this paper 
were computed on a 3D Cartesian grid of 300 3 numer- 
ical grid zones, which is about a factor three larger than 
the largest production relativity calculations of which we 
are aware (which were about 200 3 zones) . With our new 
code, these take about 12 Gbytes of memory, and require 
about a day on a 128 processor, early access Origin 2000 
supercomputer at NCSA. These calculations were carried 
out in the summer of 1997 and winter of 1998. 

B. Metric functions 

In this section we make direct comparison between the 
quantities actually being evolved in both the 2D and 3D 
codes: metric functions. This provides a direct test of 
the evolutions of the same initial data sets in completely 
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different codes, coordinate systems, and in the 3D case, 
without enforcing symmetries required in the 2D code. 



1. Geodesic slicing 



y /ip A obtained with Geodesic Slicing 



We start by considering evolutions of axisymmetric 
initial data sets with geodesic slicing (a = 1). For 
Schwarzschild, it is well known that in geodesic slicing 
a point initially on the throat will fall into the singular- 
ity at a time t = irM, and hence does not allow long 
term evolutions. For distorted black holes such as these, 
the results will differ somewhat from this analytic result, 
but the basic picture is the same. As shown in Ref. 
geodesic slicing provides a powerful tests of the evolution 
equations. For these distorted black holes, it allows us 
to compare evolved metric functions with those obtained 
from a 2D code without the complication of a lapse cal- 
culation. 

In Fig. |i| we show a plot of the conformal metric func- 
tion 7 rr = jrr/ip 4 obtained from 2D and 3D codes along 
the z-axis for the distorted black hole initial data set 
(a,b,w,n,c) = (0.5,0,1,2,0). Although the Cartesian 
metric functions are the ones actually evolved, we recon- 
struct the spherical metric functions numerically. Fur- 
thermore, dividing out the time independent conformal 
factor shows the dynamics more clearly. It is important 
to stress that the 2D data were obtained from an axisym- 
metric code that uses the logarithmic radial coordinate 
77. As discussed in section (IIIAl), this has the effect 



of providing very high resolution near the throat, unlike 
in the 3D case. The 3D data were obtained using 70 3 
grid points and a resolution of Ax = 0.0543Ma£>m- In 
Fig. we show a comparison of "fgg /r 2 from the two codes 
for the same runs. In both plots we see that the metric 
functions line up reasonably well. The small differences 
can easily be explained by the better resolution that is 
present in the 2D code, as discussed above. 

The data in Figs. [| and || are only shown to a time 
of t = Q.977Madm- The reason for this is not that the 
throat observers have already fallen into the singularity 
at that time. In fact, the 3D code with that resolution 
can be run past t — 1.95Madm without developing prob- 
lems. The problem is that the 2D code uses the standard 
spherical polar coordinate system, which has a coordinate 
singularity on the axis. Although there are no grid points 
actually on the axis, the code becomes unstable for grid 
points near the axis, causing the 2D code to crash very 
early. This problem is usually solved by using the gauge 
freedom of the Einstein equations to use the shift vector 
to maintain the metric function j r g = 0, thus keeping 
constant-coordinate lines perpendicular. However, the 
use of a non-zero shift results in a different 3-metric than 
in the vanishing shift case. Because this shift condition 
is not convenient to implement in the 3D Cartesian code, 
we were forced to compare with the zero-shift case, where 




FIG. 4. We show a comparison of ■y rr /ip 4 obtained with 
2D and 3D codes using geodesic slicing. The data are shown 
every 0.163Madm , to a maximum time of t — 0.977M adm ■ 
The 3D data were obtained with 70 3 grid points with a res- 
olution of Ax — 0.0543Madm ■ The data set evolved was 
(a,b,w,n) = (0.5,0, 1,2) (Problem 1). 



j /(r 2 ^ 4 ) obtained with Geodesic Slicing 




z/M 



FIG. 5. As in Fig. |, except that we plot jee/(r 2 tp i )- The 
data are shown at the same times as in Fig. kl 
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the 2D code develops an axis instability when the evolu- 
tion becomes very dynamic there. 

To summarize the results of this section, we have care- 
fully compared the evolutions of black hole initial data, 
in geodesic slicing, with both 2D and 3D codes. This 
provides a good test of the behavior of the code in the 
strong field regions near the horizon and singularity. In 
spite of the difficulties of evolution in 3D, in this case the 
3D evolution is more robust than the 2D case due to an 
axis instability. 



2. Singularity avoiding slicings 

We next consider evolutions with the "1+log" slic- 
ing condition. The initial lapse is obtained as described 
above in section III A 2 . Here we only show that the met- 
ric functions obtained are qualitatively the same in the 
two different codes. We cannot compare directly with 
data from 2D simulations, again because of the different 
slicing and shift conditions used. 

As an example, we consider the evolution of the initial 
data set (a,b,w,n,c) = (0.5,0,1,2,0) (Problem 1). This 
data set has an ADM mass of M ADM = 0.919M. The 
run was done with 150 3 grid points with a resolution of 
Ax — 0.054AM adm- In Fig- ^a, we show the metric func- 
tion 7 rr in the x — plane at the time t — 27.2Madm- 
We see the characteristic peak growing in this metric 
function, caused by the differential acceleration of grid 
points towards the hole. The shape of the metric func- 
tion is similar to that obtained with the 2D code. As in 
the Schwarzschild case, the peak growing near the origin 
is a result of the application of the isometry condition. 
In Fig ^|b, we show the lapse function in the x = plane 
at the same time. We see that, as in the Schwarzschild 
and Misner cases, the lapse which was initially antisym- 
metric at the throat, and evolved with "1+log" slicing, 
has collapsed at and around the throat. 

The metric functions of the other black hole spacetimes 
studied in this paper behave similarly, and the compar- 
isons with 2D results reveal the same qualitative agree- 
ment as expected. 



40 F (°) yjf 4 



C. Horizons 

In this section we show a comparison of the location 
of the apparent horizon (AH) , computed in the 3D and 
2D codes. We computed the location of the AH dur- 
ing an evolution of the initial data set (a, b, w, n, c) = 
(—0.5, 0, 1, 2, 0). The calculation was performed with 66 3 
grid points with a resolution of Ax — 0.2 = 0.056Madm- 
For this run, maximal slicing was used. The AH was 
found using the 3D AH finder detailed in Ref. E3 . This 




(b) a 




FIG. 6. We show (a) the conformal metric function 7 rr and 
(b) the lapse function a in the x = plane for the evolution 
of the initial data set (a,b,w,n,c) — (0.5,0,1,2,0) (Prob- 
lem 1). The data shown are at a time of t — 27 .2Madm ■ 
The evolution was run on a 150 3 grid, although data from 
only the inner 100 3 point are shown. The resolution was 
Ax = 0.0544M ADM . 



AH finder is implemented in full 3D, based on an expan- 
sion of the 2D horizon surface in symmetric trace free 
tensors. 

Fig. [?] shows the horizon shapes and locations in the 
3D calculation every 1.15Madm i n time, starting at t = 
0. The multipole order used in this calculation is £ = 
4, with only axisymmetric terms considered to reduce 
the computational time. The higher order expansion is 
needed here to describe the oblate shape of the horizons. 
Also shown in Fig. |?] are the corresponding surfaces found 
in our 2D axisymmetric code. We note that the slicing 
and shift conditions differ in the two cases, so we do not 
expect the surfaces to coincide precisely. In and around 
the x — y plane, the solutions match to within half a grid 
cell. Greater differences, however, are found along the z- 
axis, where the two surfaces are displaced by a maximum 
of roughly two grid cells. This is attributed in part to a 
bigger Aq (~ 0.1 ) near the z-axis, which results from 
the imposed asymmetry in the lapse function due to the 
nearness of the outer boundaries, where we enforce the 
spherical Schwarzschild lapse as a boundary condition in 
the maximal equation used in the 2D simulation. 

A more geometrical comparison or test of the solver 
is the mass of the surface found. The horizon mass 
is defined as Mah — y/ Aah/IGk, where A ah is the 
area of the surface. Fig. || plots the AH mass as a 
function of time for the 2D and the I = 4, 3D evo- 
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FIG. 7. Coordinate location of the AH found 

in the 2D and 3D evolutions of the initial data set 
(a, b, w, n, c) = (—0.5, 0, 1, 2, 0). The surfaces are shown start- 
ing at t = with time intervals of 1.12Madm ■ Although we 
do not expect identical results due to different kinematic con- 
ditions in the two cases, the surfaces differ at most by little 
more than two grid cells. The evolution is performed on a 66 3 
grid with Ax = 0.2. 

lutions. In both cases, the mass increases at first, as 
the gravitational waves fall into the black hole, reaching 
Mah ~ 0.997 M adm at t ~ QMadm- The masses in the 
two cases differ by only 0.1% at t ~ QMadm- For com- 
parison we also plot Mah for the surface found using a 
lower order 1 = 1 multipole expansion. At early times, 
when the horizon is most distorted, the t — 2 expan- 
sion is clearly not adequate to resolve the horizon shape, 
as evidenced by the AH mass which exceeds the ADM 
mass by about 1%. However, as the black hole settles 
down into a quasi-static state, the surface becomes more 
spherical and the I = 2 solution approaches both the 2D 
and the £ = 4, 3D results, differing from the 2D result by 
about 0.35% at t ~ 6Madm- 

D. Radiation extraction 

Although in black hole simulations we evolve directly 
the metric and extrinsic curvature, for applications to 
gravitational wave astronomy we are particularly inter- 
ested in computing the waveforms emitted. One mea- 
sure of this radiation is the Zerilli function, ip, which 
is a gauge-invariant function that obeys the Zerilli wave 
equation Eg ], The Zerilli function can be computed by 
writing the metric as the sum of a spherically symmet- 

o 

ric part and a perturbation: g a p —g a p +h a p, where the 
perturbation h a p is expanded in tensor spherical har- 
monics. To compute the elements of h a fj in a numerical 
simulation, one integrates the numerically evolved met- 
ric components g a p against appropriate spherical har- 
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FIG. 8. Comparison of the apparent horizon masses com- 
puted from the 2D and 3D evolutions of the initial data 
set (a,b,w,n,c) — (—0.5,0,1,2,0). The AH mass increases 
initially as the Brill wave falls into the black hole. By 
t ~ QMadm, the mass approaches 0.997Madm and the 2D 
and I = 4, 3D results differ by just 0.1% at this time. We 
also show the corresponding masses computed from the sur- 
faces found with an £ = 2 multipole expansion. The low order 
expansion is clearly not adequate in resolving the surface at 
early times when the horizon is most distorted. 

monies over a coordinate 2-sphere surrounding the black 
hole. The resulting functions can then be combined in a 
gauge-invariant way, following the prescription given by 
Moncrief j46|. This procedure was originally developed 
by Abrahams p7| , and was applied to the same class of 
distorted black hole initial data sets discussed here, but 
evolved in 2D spherical-polar coordinates and with a dif- 
ferent gauge, as discussed in ||. 

1. 3D numerical implementation 

We have developed numerical methods based on the 
same ideas to extract the waves in a full 3D Carte- 
sian setting. The method used is essentially that used 
in the axisymmctric case, except that the metric func- 
tions and their spatial derivatives need to be interpo- 
lated onto a two-dimensional surface, which we choose to 
have constant coordinate radius. The projections of the 
perturbed metric functions h a p, and their radial deriva- 
tives, are then computed by numerically performing two- 
dimensional surface integrals for each £ — m mode de- 
sired. Then, for each mode, the Zerilli function is con- 
structed from these projected metric functions, accord- 
ing to Moncrief 's gauge-invariant prescription. This is a 
complicated but straightforward procedure. 

a. Constructing the Zerilli Function As mentioned 
above, we assume the general metric can be decomposed 
into its spherical and non-spherical parts: 
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(8) 



The spherical part g^ will of course be Schwarzschild, 
but we will in general not know the mass of this 
Schwarzschild background, or what coordinate system it 
will be in. However, in general, we know it can be written 
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where the functions N, A, and R are functions of our 
coordinate radius and time. Regge and Wheeler showed 
that hfu, for even-parity perturbations (i.e., perturba- 
tions which do not introduce angular momentum) can be 
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The spherical part of the metric, given in the functions 
N, A, and R, can be obtained by projecting the full met- 
ric against Y 00 , yielding the following expressions^]: 
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(20) 
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Each £m-mode of can then be obtained by projecting 
the full metric against the appropriate Yg m . For example: 



H, 
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(23) 



*We thank Gabrielle Allen for pointing out the error in the 
expression for R 2 in reference M. 



Expressions for the other functions are provided in 
Ref. |22j] . In practice, we do not extract A to compute 
V°ut rather we assume it to have the form 1 — ^ , 
where we take M to be the ADM mass of the spacetime. 
For the axisymmetric, equatorial plane symmetric, non- 
rotating initial data studied here, the only modes allowed 
are even-parity, even-£ modes. 

Moncrief showed that the Zcrilli function is gauge in- 
variant, and can be constructed from the Regge- Wheeler 
variables as follows: 



On) 



where 



h{£- 1)0? + 2) 4RS 2 k ( 2 lm) 



-£{£+l)Rk\ 



{em) 



£(£ + !) 



A 



A 



£(£- 



SRG 



6M 
~R 

{im 



R 



-2-h {l 
R hl 



lira) 



H. 



{im) 



2S 



1 d (RK (ln 
2~7^dR \ 



S = l 



2M 



(24) 

(25) 
(26) 

(27) 

(28) 



In order to compute the Regge- Wheeler perturbation 
functions hi, H2, G, and K, one needs the spherical met- 
ric functions on some 2-sphere. We get these by inter- 
polating the Cartesian metric functions onto a surface of 
constant coordinate radius, and computing the spherical 
metric functions from these using the standard transfor- 
mation. Second order interpolation is used. In order to 
compute the needed radial derivatives of these functions, 
we first compute the derivative of the Cartesian metric 
functions with respect to the coordinate radius on the 
Cartesian grid. We interpolate these quantities onto the 
2-sphere. Then, to get the derivatives with respect to the 
Schwarzschild radius R, we use the derivative of Eq. (I2T 
with respect to the coordinate r, giving 



dR 
dr 



1 



IQttR 



gee,r 



gc/>c/>,r 

sin 2 9 



dn. 



(29) 



The points on the 2-sphere to which we interpolate 
lie on lines of constant 9 and </>, staggered across the 
coordinate axes. Integrations over each of these variables 
is performed with the second-order accurate Simpson's 
rule. Currently, 300 points are used for each angular 
coordinate. The numerical interpolations involved in this 
extraction procedure were also chosen to be second order 
accurate. Testing shows that both the interpolations and 
integrations do converge to second order in the relevant 
grid spacing. 

We are aware of two possible improvements on this 
scheme. First, in the ^-integration, we could integrate 
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over cos 9 instead of 9. This could make these integra- 
tions more accurate, and would allow us to use fewer 
points, making the code more efficient. Second, we could 
interpolate onto a stereographic patch. This would dis- 
tribute the points more evenly on the surface. 

As in Ref. || , we choose to normalize the Zerilli func- 
tion so that the asymptotic energy flux in each mode is 
given by 

E = (1/32tt)V; 2 . (30) 

This assumes that linear theory is adequate to treat each 
individual mode. However, it is possible that in some 
regimes linear theory is not adequate, and hence the 
true energy calculation would require additional terms. 
The treatment of second order perturbation theory of 
this sort, including a derivation of the energy formula in- 
cluding higher order corrections, has been developed in 
Ref. @. 

Such considerations are actually relevant for the some 
of the initial data sets described here. For example, Prob- 
lem 1 has angular parameter n = 2. One can show |3C| l 
that in an expansion of the initial data for small ampli- 
tude a, one recovers the Schwarzschild background with 
nonspherical perturbations, as assumed above. However, 
for n = 2, the series expansion in the amplitude parame- 
ter a brings in the £ = 2 perturbations at linear order, but 
£ = 4 does not appear until second order in a. Therefore, 
in order to be fully self-consistent the evolution equa- 
tions for the £ — 4 modes should actually be treated to 
second order, which implies that there should be nonlin- 
ear source terms coming from other linear modes (in this 
case, £ = 2) in a consistent perturbation expansion, as 
discussed in Ref. Q. Basically, in such cases the £ = 4 
mode is so small, that nonlinear contributions from other 
modes are not to be neglected. This need not concern 
us, because the evolutions considered in this paper are 
carried out fully nonlinearly, but it does mean that the 
interpretation, and in particular the calculation of the 
energy carried by a mode, must be carefully considered. 
Perturbative calculations, with comparisons to full non- 
linear simulations, are in progress and will be reported 
elsewhere ||^3|. 

In the results discussed below, in cases where second 
order corrections should be considered in principle, we 
will continue to use the linear energy formula as a non- 
rigorous indication of the strength of a signal. The cases 
in question concern the £ > n modes (e.g., £ — 6 in all 
cases and £ = 4 for the n = 2 simulations, where n is the 
angular index of the Brill wave.) 

To summarize this section, we have developed a gen- 
eral technique to extract the 3D waves propagating on a 
black hole background. There are cases where linear the- 
ory must be used with caution. While previously only 
axisymmetric simulations have been studied, we can now 
study all non-trivial wave modes, including those with 
to 7^ 0. 



2. Waveforms from 3D simulations of distorted black holes 

In this section, we examine the waveforms for the three 
test cases described in detail above. This is a detailed 
followup to Ref. |]lj, so we present £ — 2 and £ = 4 
waveforms for Problem 1 as before. We also consider two 
additional cases with different amplitudes and angular 
index n to show the fairly generic nature of these results. 
Finally, for the first time, we show that it is now possible 
in some cases to extract accurately even the £ = 6 even- 
parity waveforms, both in the 2D and 3D cases, if the 
amplitude is not too small. 

a. Problem 1. We extracted the £ = 2 and £ = 4 
Zerilli functions during an evolution of the distorted black 
hole initial data set (a,b,w,n,c) — (0.5,0,1,2,0), using 
the extraction method described above. In Fig. ||a we 
show the £ = 2 Zerilli function extracted at a radius r = 
8.16Madm as a function of time. Superimposed on this 
plot is the same function computed during the evolution 
of the same initial data set with a 2D code, based on 
the one described in detail in ||J3^]. The agreement of 
the two plots over the first peak is a strong affirmation 
of the 3D evolution code and extraction routine. It is 
important to note that the 2D results were computed 
with a different slicing (maximal), different coordinate 
system, and a different spatial gauge. Yet the physical 
results obtained by these two different numerical codes, 
as measured by the waveforms, are remarkably similar (as 
one would hope). This is a principal result of this paper. 
A full evolution with the 2D code to t — IOOMadm, by 
which time the hole has settled down to Schwarzschild, 
shows that the energy emitted in this mode at that time 
is about 4 x 10~ 3 Madm- Therefore, although this was a 
highly distorted black hole, less than 5% of the MRL is 
actually radiated away. (Other modes are much smaller 
and do not contribute significantly to the total energy 
radiated.) 

In Fig. |^b we show the £ — 4 Zerilli function extracted 
at the same radius, computed during evolutions with 2D 
and 3D codes. This waveform is more difficult to extract, 
because it has a higher frequency in both its angular 
and radial dependence, and it has a much lower ampli- 
tude: using the linear theory formula discussed above 
(Eq. (j30|)) the energy emitted in this mode is three or- 
ders of magnitude smaller than the energy emitted in 
the £ = 2 mode, i.e., 10~ 6 Madm, yet it can still be ac- 
curately evolved and extracted. 

Small differences between the 2D and 3D results can 
be seen. Resolution studies of the 3D results indicate 
that the differences are not completely due to resolution 
of the 3D evolution code. The small differences in phase 
can be understood as a result of the different shift and 
slicings being used in the two simulations. The radia- 
tion is extracted at a constant coordinate location, and 
the coordinates fall towards the black hole at different 
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FIG. 9. We show the (a) 1 = 2, and (b) I = 4 Zerilli func- 
tions vs. time, extracted during 2D and 3D evolutions of the 
data set (a, b, w, n, c) = (0.5, 0, 1, 2, 0) (Problem 1). The func- 
tions were extracted at a radius of 8.19Madm ■ The 2D data 
were obtained with 202 x 54 grid points, giving a resolution 
of Arj = A6 = 0.03. The 3D data were obtained using 300 3 
grid points and a resolution of Ax = 0.109Madm ■ 



rates with different slicings and shifts. By measuring the 
physical radial position of the wave extraction in these 
simulations, we determined that the difference between 
the 2D and 3D phases at late time is consistent with the 
slightly different extraction locations in the two cases. 
The additional differences in the I — 4 waveforms could 
be related to slight differences in the initial data, which 
were generated in independent ways, or even differences 
in gauge (the waveforms are gauge-invariant, meaning 
they are unaffected only at first order under gauge trans- 
formations). As £ — 4 has a much smaller amplitude 
than £ = 2, it will be more sensitive to such details. 

Since this data set is axisymmetric, the non- 
axisymmctric Zerilli functions of a perfect evolution 
should vanish. However, we expect that numerical errors, 
especially errors due to the distinctly non-axisymmetric 
outer boundary, will result in finite values for these Zer- 
illi functions. Knowing the size of the error produced 
will help us gauge the accuracy of our results when we 
move to full 3D, where we expect real non-axisymmetric 
signals. 

In Fig. [l^, we show non-axisymmetric Zerilli functions 
as a function of time, extracted during an evolution of 
the same data set as above, at the same radius. We see 
that the £ = 2, m = 2 and £ = 4, m = 2 waveforms 
remain very small throughout the evolution. The £ — 4, 
m = 4 waveform, while an order of magnitude smaller 
than the £ = 4, m = waveform, is still much larger 
than the other non-axisymmetric waveforms. Prelimi- 
nary tests computing these quantities for Schwarzschild 
initial data, for which all Zerilli functions should vanish, 
indicate that the £ = 4, m = 4 waveform is more sen- 
sitive to errors in the interpolation of metric functions 
onto the spherical extraction surface than the other Zer- 
illi functions are. The reason for this sensitivity and how 
to correct for it are under investigation. 

Resolution studies show that the waveforms are con- 
verging at roughly second order in the grid spacing. The 
axisymmetric £ — 2, 4 modes have converged at the res- 
olution shown, while the £ = m = 4 mode is converging 
towards zero, although at the highest resolution we have 
computed it is still relatively large. 

b. Problem 2 We now turn to Problem 2, which is 
similar to Problem 1, but with angular index n = 4. In 
this case the initial data contain a much larger £ — 4 
and £ = 6 mode content. In this example, we extract the 
£ = 2,4,6 modes during the 3D evolution, and compare 
with the 2D axisymmetric evolution. 



Fig. 
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shows the 



comparisons between the waveforms, showing that for 
this highly distorted black hole, even the £ = 6 waveform 
computed in 3D agrees very well with the axisymmetric 
simulation. As in the previous case, we emphasize that 
these two calculations are done in different geometries, 
different slicing conditions, and different spatial gauges. 

For this simulation, the total radiated energy carried 
in the £ = 2 mode is 2 x W~ 3 Madm , or about 3% of the 
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FIG. 10. We show the non-axisymmetric Zerilli functions 
as a function of time extracted during an evolution of the 
axisymmetric data set (a, b, w, n, c) = (0.5, 0, 1, 2, 0) (Problem 
1). The extraction was done at a radius of r — 8.16Madm ■ 
The calculation was performed with 300 3 grid points and a 
resolution of Ax = 0.109Madm ■ 

MRL. The £ = 4 mode carries 1 x 10~ 4 Madm in energy, 
much larger than in the previous case. In this case the 
£ = 4 mode is a first order perturbation in the Brill wave 
amplitude, and the linear energy calculation applies. 

c. Problem 3 Finally, we study the waveforms ex- 
tracted in Problem 3, which has a much lower distor- 
tion amplitude, and hence much weaker waves which are 
harder to evolve and extract accurately. In Fig. [jj] we 
show the three waveforms. The £ = 2 and £ = 4 3D 
waveforms agree well with those computed in 2D, as in 
previous cases. But the £ = 6 waveform has an offset 
indicating some level of error in cither the extraction or 
the evolution itself. In this simulation the £ = 6 mode is 
about 20 times smaller than that extracted in Problem 
2. If one examines the £ = 6 waveform from Problem 
2 closely, one can see a similar offset of the same level, 
but compared to the waveform in that case the effect 
is almost not noticeable. However, in Problem 3, with 
such a smaller amplitude signal, the offset can be seed 
readily. Resolution studies of this mode indicate that it 
is converging toward the 2D result as the resolution in- 
creases. But even at the highest 300 3 resolution, there 
is enough error to make extracting this very small signal 
(much smaller than any other modes presented) rather 
difficult. The reasons for the offset in £ = 6 waveforms 
are under investigation. 

The energy emitted in the £ = 2 mode is 9 x 
10~ 5 Madm, again, only a few per cent of the MRL. The 
energy emitted in the £ = 4 mode is 6 x 10~ s Madm- 



IV. CONCLUSIONS 
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FIG. 11. We show the (a) i = 2, (b) I = 4, and (c) I = 6 
Zerilli functions vs. time, extracted during 2D and 3D evolu- 
tions of the data set (a,b,w,n,c) = (0.5,0,1,4,0) (Problem 
2). The functions were extracted at a radius of 7.72Madm- 
The 2D data were obtained with 202 x 54 grid points, giving 
a resolution of A-q = A8 = 0.03. The 3D data were obtained 
using 300 3 grid points and a resolution of Ax — 0.103Madm- 



In this followup paper to Ref. JlJ, we have performed 
the first detailed study of the dynamics of axisymmetric, 
distorted black holes with a general 3D code. This work 
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FIG. 12. We show the (a) I = 2, (b) I = 4, and (c) ^ = 6 
Zerilli functions vs. time, extracted during 2D and 3D evolu- 
tions of the data set (a,b,w,n,c) = (0.1,0, 1,4,0) (Problem 
3). The functions were extracted at a radius of 7.88Madm- 
The 2D data were obtained with 202 x 54 grid points, giving 
a resolution of A-q = A9 — 0.03. The 3D data were obtained 
using 300 3 grid points and a resolution of Ax = O.IOAMadm- 



is another step in the development of 3D numerical rel- 
ativity towards codes capable of simulating the spiraling 
coalescence of binary black hole systems. In this step, we 
have shown that, given sufficient resolution in 3D Carte- 
sian coordinates, distorted black holes, which are used to 
model the late stages of binary black hole mergers just af- 
ter the horizons have merged, can be accurately evolved. 
We studied the behavior of metric functions and apparent 
horizons in these spacetimes evolved in 3D Cartesian co- 
ordinates, and compared them with evolutions provided 
by 2D codes in polar-spherical coordinates, and showed 
that the results agree well. 

Furthermore, we showed that the gravitational wave- 
forms generated by the black hole, consisting of small 
perturbations on the evolving black hole background, can 
be accurately propagated and extracted from the numeri- 
cally generated metric, on a 3D Cartesian grid. We have 
demonstrated this by comparing results from a mature 
2D code, showing good agreement not only for the I = 2, 
but also with the higher I — 4 and even the t = 6 modes 
of the radiation, in cases where they are strong enough 
to stand above numerical error. Such modes are much 
weaker than the dominant £ = 2 modes, and have much 
higher frequency and more complex angular dependence, 
yet they can be accurately tracked and resolved in full 
3D. To our knowledge, this is the first study of £ — 6 
modes with any nonlinear evolution code. 

Although we regard this as an important step in es- 
tablishing numerical relativity as a viable tool to com- 
pute waveforms from black hole interactions, the calcu- 
lations one can presently do are limited. With present 
techniques, the evolutions can only be carried out for a 
fraction of the time required to simulate the 3D orbit- 
ing coalescence. Many techniques to handle this more 
general case are under development, such as hyperbolic 
formulations of the Einstein equations and the advanced 
numerical methods they bring |fl9f , adaptive mesh refine- 
ment that will enable placing the outer boundary farther 
away while resolving the strong field region where the 
waves are generated | p0| , and apparent horizon bound- 
ary conditions that excise the interiors of the black holes, 
thus avoiding the difficulties associated with singularity 
avoiding slicings @|||jl|] . 

In a completely different approach, a 3D code based 
on a characteristic formulation of the equations was 
shown to overcome some limitations of 3D Cauchy codes, 
to evolve distorted black holes for essentially unlimited 
times (more than t « 60, 000M), on distorted black holes. 
Although it is not clear if characteristic evolution will be 
able to handle highly distorted or colliding black holes, 
it would be very interesting to study black hole systems 
such as these to see if the black hole ringdown and wave- 
forms can be also be accurately extracted. 

All of these techniques, and others, may be needed 
to handle the more general, long term evolution of coa- 
lescing black holes. Our purpose in this paper has been 
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to show that (a) given present resources one can evolve 
simpler distorted black hole systems and accurately ex- 
tract the waveforms, even when they carry only 10 _6 M 
in energy, and (b ) to establish testbeds for the techniques 
under development for the more general case. Experience 
has shown that although the gross dynamics of black hole 
evolutions are fairly robust to different treatments in fully 
nonlinear codes, the waveforms are very sensitive and dif- 
ficult to compute correctly |Q,0-^2|j2|Q . Each of these 
new techniques may introduce numerical artifacts, even 
if at very low amplitude, to which the waveforms may be 
very sensitive. As new methods are developed and ap- 
plied to numerical black hole simulations, they can now 
be tested on evolutions such as those presented here to 
ensure that the waveforms can be accurately computed. 

In future papers we extend this work and apply it to 
full 3D initial data sets where nonaxisymmetric modes 
can be extracted for the first time pi] , p3[j52|| , and to the 
evolution of colliding black holes in 3D, extending the 
work in |io). Once these techniques have been fully de- 
veloped and tested on true 3D data sets, it will be im- 
portant to apply it to true 3D black hole collision sim- 
ulations, such as those recently reported by Brugmann 
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